home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / tri_surf.pro < prev    next >
Text File  |  1997-07-08  |  7KB  |  201 lines
  1. ; $Id: tri_surf.pro,v 1.5 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1993-1997, Research Systems, Inc.  All rights reserved.
  4. ;    Unauthorized reproduction prohibited.
  5.  
  6. FUNCTION TRI_SURF, z, x, y, REGULAR = regular, XGRID=xgrid, $
  7.     XVALUES = xvalues, YGRID = ygrid, YVALUES = yvalues, $
  8.     GS = gs, BOUNDS = bounds, NX = nx0, NY = ny0, $
  9.     MISSING = missing, EXTRAPOLATE = extrapolate, LINEAR=linear
  10. ;+
  11. ; NAME:
  12. ;    TRI_SURF
  13. ;
  14. ; PURPOSE:
  15. ;    This function Interpolates a regularly or irregularly gridded
  16. ;    set of points with a smooth quintic surface.
  17. ;
  18. ; CATEGORY:
  19. ;    Mathematical functions, Interpolation, Surface Fitting
  20. ;
  21. ; CALLING SEQUENCE:
  22. ;    Result = TRI_SURF(Z [, X, Y])
  23. ;
  24. ; INPUTS: 
  25. ;    X, Y, Z:  arrays containing the X, Y, and Z coordinates of the 
  26. ;          data points on the surface. Points need not be 
  27. ;          regularly gridded. For regularly gridded input data, 
  28. ;          X and Y are not used: the grid spacing is specified 
  29. ;          via the XGRID and YGRID (or XVALUES and YVALUES) 
  30. ;          keywords, and Z must be a two dimensional array. 
  31. ;          For irregular grids, all three parameters must be
  32. ;          present and have the same number of elements. 
  33. ;
  34. ; KEYWORD PARAMETERS:
  35. ;  Input grid description:
  36. ;    REGULAR:  if set, the Z parameter is a two dimensional array 
  37. ;          of dimensions [N,M], containing measurements over a 
  38. ;          regular grid. If any of XGRID, YGRID, XVALUES, YVALUES 
  39. ;          are specified, REGULAR is implied. REGULAR is also 
  40. ;          implied if there is only one parameter, Z. If REGULAR is 
  41. ;          set, and no grid (_VALUE or _GRID) specifications are 
  42. ;          present, the respective grid is set to (0, 1, 2, ...). 
  43. ;    XGRID:    contains a two element array, [xstart, xspacing], 
  44. ;          defining the input grid in the X direction. Do not
  45. ;          specify both XGRID and XVALUES. 
  46. ;    XVALUES:  if present, XVALUES[i] contains the X location 
  47. ;          of Z[i,j]. XVALUES must be dimensioned with N elements. 
  48. ;    YGRID:    contains a two element array, [ystart, yspacing], 
  49. ;          defining the input grid in the Y direction. Do not
  50. ;          specify both YGRID and YVALUES. 
  51. ;    YVALUES:  if present, YVALUES[i] contains the Y location 
  52. ;          of Z[i,j]. YVALUES must be dimensioned with N elements. 
  53. ;
  54. ;  Output grid description:
  55. ;    GS:      If present, GS must be a two-element vector [XS, YS],
  56. ;          where XS is the horizontal spacing between grid points
  57. ;          and YS is the vertical spacing. The default is based on
  58. ;          the extents of X and Y. If the grid starts at X value
  59. ;          Xmin and ends at Xmax, then the default horizontal
  60. ;          spacing is (Xmax - Xmin)/(NX-1). YS is computed in the
  61. ;          same way. The default grid size, if neither NX or NY
  62. ;          are specified, is 26 by 26. 
  63. ;    BOUNDS:   If present, BOUNDS must be a four element array containing
  64. ;          the grid limits in X and Y of the output grid:
  65. ;          [Xmin, Ymin, Xmax, Ymax]. If not specified, the grid
  66. ;          limits are set to the extent of X and Y. 
  67. ;    NX:       The output grid size in the X direction. NX need not
  68. ;            be specified if the size can be inferred from GS and
  69. ;          BOUNDS. The default value is 26.
  70. ;    NY:       The output grid size in the Y direction. See NX. 
  71. ;
  72. ;  Others:
  73. ;       EXTRAPOLATE: If set, extrapolate the surface to points outside
  74. ;        the convex hull of input points.  Has no effect if
  75. ;        input points are regularly gridded.
  76. ;    MISSING: If set, points outside the convex hull of input points
  77. ;        are set to this value.  Default is 0.  Has no effect
  78. ;        if input points are regularly gridded.
  79. ;    LINEAR:  If set, linear interpolation, rather than quintic is
  80. ;        used.  Gradient estimates are not used.  This is faster
  81. ;        and numerically stable, but does not extrapolate.
  82. ;
  83. ; OUTPUTS:
  84. ;    This function returns a two-dimensional floating point array
  85. ;     containing the interpolated surface, sampled at the grid points.
  86. ;
  87. ; RESTRICTIONS:
  88. ;    The output grid must enclose convex hull of the input points.
  89. ; PROCEDURE:
  90. ;    This routine is similar to MIN_CURVE_SURF but the surface
  91. ;    fitted is a smooth surface, not a minimum curvature surface.  This
  92. ;    routine has the advantage of being much more efficient
  93. ;    for larger numbers of points.
  94. ;
  95. ;    The built-in IDL routines TRIANGULATE and TRIGRID(/QUINTIC) are
  96. ;    used.
  97. ;
  98. ; EXAMPLES:
  99. ; Example 1: Irregularly gridded cases
  100. ;    Make a random set of points that lie on a gaussian:
  101. ;      n = 15                ;# random points
  102. ;      x = RANDOMU(seed, n)
  103. ;      y = RANDOMU(seed, n)
  104. ;      z = exp(-2 * ((x-.5)^2 + (y-.5)^2))      ;The gaussian
  105. ;
  106. ;     get a 26 by 26 grid over the rectangle bounding x and y:
  107. ;      r = TRI_SURF(z, x, y)    ;Get the surface.
  108. ;
  109. ;     Or: get a surface over the unit square, with spacing of 0.05:
  110. ;      r = TRI_SURF(z, x, y, GS=[0.05, 0.05], BOUNDS=[0,0,1,1])
  111. ;
  112. ;     Or: get a 10 by 10 surface over the rectangle bounding x and y:
  113. ;      r = TRI_SURF(z, x, y, NX=10, NY=10)
  114. ;
  115. ; Example 2: Regularly gridded cases
  116. ;    Make some random data
  117. ;      z = randomu(seed, 5, 6)
  118. ;
  119. ;    interpolate to a 26 x 26 grid:
  120. ;      CONTOUR, TRI_SURF(z, /REGULAR)
  121. ;
  122. ; MODIFICATION HISTORY:
  123. ;    DMS, RSI, October, 1993.  Written.
  124. ;    DMS, RSI, Jan, 1996. Added LINEAR keyword.
  125. ;-
  126.  
  127.  
  128. ON_ERROR, 2
  129. s = size(z)        ;Assume 2D
  130. nx = s[1]
  131. ny = s[2]
  132.  
  133. reg = keyword_set(regular) or (n_params() eq 1)
  134.  
  135. if n_elements(xgrid) eq 2 then begin
  136.     x = findgen(nx) * xgrid[1] + xgrid[0]
  137.     reg = 1
  138. endif else if n_elements(xvalues) gt 0 then begin
  139.     if n_elements(xvalues) ne nx then $
  140.         message,'Xvalues must have '+string(nx)+' elements.'
  141.     x = xvalues
  142.     reg = 1
  143. endif
  144.  
  145. if n_elements(ygrid) eq 2 then begin
  146.     y = findgen(ny) * ygrid[1] + ygrid[0]
  147.     reg = 1
  148. endif else if n_elements(yvalues) gt 0 then begin
  149.     if n_elements(yvalues) ne ny then $
  150.         message,'Yvalues must have '+string(ny)+' elements.'
  151.     y = yvalues
  152.     reg = 1
  153. endif
  154.  
  155. if reg then begin
  156.     if s[0] ne 2 then message,'Z array must be 2D for regular grids'
  157.     if n_elements(x) ne nx then x = findgen(nx)
  158.     if n_elements(y) ne ny then y = findgen(ny)
  159.     x = x # replicate(1., ny)    ;Expand to full arrays.
  160.     y = replicate(1.,nx) # y
  161.     endif
  162.  
  163. n = n_elements(x)
  164. if n ne n_elements(y) or n ne n_elements(z) then $
  165.     message,'x, y, and z must have same number of elements.'
  166.  
  167. if n_elements(bounds) lt 4 then begin    ;Bounds specified?
  168.     xmin = min(x, max = xmax)
  169.     ymin = min(y, max = ymax)
  170.     bounds = [xmin, ymin, xmax, ymax]
  171.     endif
  172.  
  173. TRIANGULATE, x, y, tr, b
  174.  
  175. if n_elements(gs) lt 2 then begin    ;GS specified?  No.
  176.     if n_elements(nx0) le 0 then nx = 26 else nx = nx0 ;Defaults for nx and ny
  177.     if n_elements(ny0) le 0 then ny = 26 else ny = ny0
  178.     gs = [(bounds[2]-bounds[0])/(nx-1.), $
  179.        (bounds[3]-bounds[1])/(ny-1.)]
  180.   endif else begin            ;GS is specified?
  181.     if n_elements(nx0) le 0 then $
  182.     nx = ceil((bounds[2]-bounds[0])/gs[0]) + 1 $
  183.     else nx = nx0
  184.     if n_elements(ny0) le 0 then $
  185.     ny = ceil((bounds[3]-bounds[1])/gs[1]) + 1 $
  186.     else ny = ny0
  187.   endelse
  188.  
  189.  
  190. if N_ELEMENTS(missing) le 0 then missing = 0
  191.  
  192. if KEYWORD_SET(linear) then $
  193.     return, TRIGRID(x,y,z, tr, gs, bounds, MISSING=missing)
  194.  
  195. if KEYWORD_SET(extrapolate) then $
  196.     return, TRIGRID(x,y,z, tr, /QUINTIC, gs, bounds, $
  197.             MISSING=missing, EXTRAPOLATE=b)
  198.  
  199. return, TRIGRID(x,y,z, tr, /QUINTIC, gs, bounds, MISSING=missing)
  200. end
  201.